/*
    Linear Programming Library
    
    Created by Team Paragon (Shena Hoffmann, Arthur Powalka, Nasir Kenarangui)
    CSCE 689
    
    RevSimplex.cpp
    Defines the variables and functions used for providing functionality related to the revised simplex linear programming algorithm
*/

#include "RevSimplex.h"

using namespace arma;


// ............ CONSTRUCTORS ............ //
RevSimplex::RevSimplex()
{ }

RevSimplex::RevSimplex(const LinProg &lp) : LinProg(lp)
{ }

RevSimplex::RevSimplex(const RevSimplex &lp) : LinProg(lp)
{ }

RevSimplex::RevSimplex(bool max, vector<double> objective, vector<vector<double>> constraints, vector<string> inequalitySigns, vector<double> inequalityVals, vector<bool> constrainedVars)
{
    this->max = max;
    this->objective = objective;
    this->constraints = constraints;
    this->inequalitySigns = inequalitySigns;
    this->inequalityVals = inequalityVals;
    this->constrainedVars = constrainedVars;
    verifyLinProg();
}

// ............ DESTRUCTOR ............ //

RevSimplex::~RevSimplex()
{
}

// ............ FUNCTIONS FOR THE USER ............ //

//NOTE: all credit for the algorithm in this function goes to Dr. Kartik at North Carolina State University (http://www4.ncsu.edu/~kksivara/ma505/codes/revised_simplex.m)
tuple<vector<double>, string, double> RevSimplex::solve()
{ 
	LinProg standardForm=convertToStandardForm();	
	RevSimplex test(standardForm);	
	RevSimplex final=test.convertToRsimplexForm();
	setClassVariables(final);
	C=trans(C);
	double eps=pow(10,-3);
	mat temp;	
	mat tempB;
	for(int k=standardForm.getObjective().size();k<A.n_cols;k++)	//will determine the initial basis. The indices correspond to the columns which correspond to the variables
	{
		temp=k;
		tempB=join_cols(tempB, temp);
	}
	B=conv_to<uvec>::from(tempB);

	uvec N(A.n_cols);
	
	//creates a vector with nonbasic variable indices
	for(int e=0;e<A.n_cols;e++)
	{
		N(e)=e;	
	}	
	N.shed_rows(B(0),B(B.n_elem-1));
	mat xB=inv(A.cols(B))*b;
	while(true)
	{
		mat cB=C.rows(B);
		mat y=inv(trans(A.cols(B)))*cB;
		colvec sN=C.rows(N)-trans(A.cols(N))*y;
		double sNMax=arma::max(sN);
		string feasible;	
		vector<double> varSolution;	
		double objValue=0;
		if (sNMax<=eps)
		{
			mat x=zeros<mat>(A.n_cols,1);
			x.rows(B)=xB;
			bool testFeasibility=1;
			for(int f=0;f<x.n_rows;f++) //creates solution vector for the variables
			{
				varSolution.push_back(x(f));
				if(x(f)<0)
				{
					testFeasibility=0; //a negative solution cannot occur
				}	
			}
			if(testFeasibility)
			{
				feasible="feasible";
				mat zValue=trans(C)*x;
				objValue=zValue(0,0);
				
			}
			else	
			{
				feasible="not feasible";
				varSolution.clear();
			}
			auto solutionData=make_tuple(varSolution, feasible, objValue);
			C.clear();
			A.clear();
			b.clear();
			B.clear();
			return solutionData;
		}
		uvec index=find(sN==sNMax); //don't use equality
		int k=index(0);
		//finding the leaving variable
		mat d=inv(A.cols(B))*A.col(N(k));
		uvec zz=find(d>eps);	
		if(zz.is_empty())
		{
			feasible="not feasible or unbounded";
			mat x=zeros<mat>(A.n_cols,1);
			x.rows(B)=xB;
			auto solutionData=make_tuple(varSolution, feasible, objValue);
			C.clear();
			A.clear();
			b.clear();
			B.clear();
			return solutionData;
		}

		mat test=xB.rows(zz)/d.rows(zz);
		mat theta=min(test);
		double theta_doub=theta(0);
		uvec ii=find(test==theta_doub);
	
		double l=zz(ii(0));
		int temp3=B(l); 
		B(l)=N(k);
		N(k)=temp3;
		xB=xB-theta_doub*d;

		xB(l)=theta_doub;
	
	}

	C.clear();
	A.clear();
	b.clear();
	B.clear();
}

template<typename type>
void RevSimplex::printVector(vector<type> vec)	
{
	for_each(vec.begin(), vec.end(), [](type value){cout<<value<<"\t";});
	cout<<endl;
}

vector< complex<double> > RevSimplex::convertObjectiveToComplexForm(vector<double> oldObjective)
{
    vector< complex<double> > newComplexObjective;
    for (int i = 0 ; i <oldObjective.size(); i++)
    {
        newComplexObjective.push_back(oldObjective[i]);
    }
    return newComplexObjective;
}

RevSimplex RevSimplex::convertToRsimplexForm()
{
    verifyLinProg();
    RevSimplex StdSmx;
    StdSmx.complexObjective = StdSmx.convertObjectiveToComplexForm(objective);    
    StdSmx.objective = objective;
    StdSmx.constraints = constraints;
    StdSmx.inequalitySigns = inequalitySigns;
    StdSmx.inequalityVals = inequalityVals;
    StdSmx.constrainedVars = constrainedVars;
    int Osize = StdSmx.complexObjective.size();
    int Csize = StdSmx.constraints.size();
    double k =0;
    complex <double> l= 0;
    //note we do not flip the signs in this version. Simplex has a version that has an extra for loop here
    // change min to max
    if (max == false)
    {
        for (int i = 0 ; i < Osize; i++)
        {
            l = StdSmx.complexObjective[i];
           StdSmx.complexObjective[i] = -l;
        }
    }
    
    // flip sign by changing inequality
    k=0;
    
    // add slack
    for(int i=0;i<Csize;i++)	//equation i
    {
        if(StdSmx.inequalitySigns[i]=="=")
        {
            StdSmx.constraints[i].push_back (1);
            StdSmx.complexObjective.push_back({0,-1});
            for(int k=0;k<Csize;k++)
            {
                if( k!= i)
                {
                    StdSmx.constraints[k].push_back (0);		
                }
            }
        }
        else if(StdSmx.inequalitySigns[i]=="<=")
        {
            StdSmx.constraints[i].push_back(1);
            StdSmx.complexObjective.push_back(0);
            for(int k=0;k<Csize;k++)
            {
                if( k!= i)
                {
                    StdSmx.constraints[k].push_back (0);		
                }  
            }
        }
        else if(StdSmx.inequalitySigns[i]==">=")
        {
            StdSmx.constraints[i].push_back(-1);
            StdSmx.constraints[i].push_back(1);
            StdSmx.complexObjective.push_back(0);
            StdSmx.complexObjective.push_back({0,-1});
            for(int k=0;k<Csize;k++)
            {
                if( k!= i)
                {
                    StdSmx.constraints[k].push_back (0);	
                    StdSmx.constraints[k].push_back (0);	
                }  
            }
        }
    }
    
    // change to sign constrained
    for (int i = 0 ; i < Osize; i++)
    {
        if (StdSmx.constrainedVars[i]!=true)
        {	
            for (int j =0; j<Csize ; j++)
            {
                if (StdSmx.constraints[j][i]!=0)
                {
                    k = StdSmx.constraints[j][i];
                    StdSmx.constraints[j].push_back(-k);
                }
                else
                    StdSmx.constraints[j].push_back(0);

            }
            StdSmx.complexObjective.push_back(-StdSmx.complexObjective[i]);
        }
    }
    k=0;
    return StdSmx;
}

void RevSimplex::setClassVariables(LinProg &obj)
{
	rowvec temp = conv_to<rowvec>::from(obj.getObjective());
	C=temp;
	vector< vector<double> > conVec=obj.getConstraints();

	for_each(conVec.begin(),conVec.end(), [&](vector<double> value)
	{
		rowvec temp = conv_to<rowvec>::from(value);
		A=join_cols(A,temp);
	});

	b = conv_to<mat>::from(obj.getInequalityVals());
}

void RevSimplex::setClassVariables(RevSimplex &obj) 
{
	vector<double> smallObjective=obj.getObjective();
	int n=smallObjective.size();
	int m=(obj.getConstraints())[0].size();
	for(int y=0;y<(m-n);y++)
	{
		smallObjective.push_back(0);
	}
	rowvec temp = conv_to<rowvec>::from(smallObjective);
	C=temp;
	vector< vector<double> > conVec=obj.getConstraints();

	for_each(conVec.begin(),conVec.end(), [&](vector<double> value)
	{
		rowvec temp = conv_to<rowvec>::from(value);
		A=join_cols(A,temp);
	});

	b = conv_to<mat>::from(obj.getInequalityVals());
}

// ............ OVERLOADED OPERATORS ............ //

RevSimplex RevSimplex::operator=(const RevSimplex &lp)
{
    max = lp.max;
    objective = lp.objective;
    constraints = lp.constraints;
    inequalitySigns = lp.inequalitySigns;
    inequalityVals = lp.inequalityVals;
    constrainedVars = lp.constrainedVars;
}
